Structural Susceptibility and Separation of Time Scales in the van der Pol Oscillator 



<N 

O 
(N 

bJQ[ 



Ricky Chachra, Mark K. Transtrum, and James P. Sethnalj 
Department of Physics, 
Laboratory of Atomic and Solid State Physics 
Cornell University, Ithaca, NY 14853 USA 
(Dated: August 3, 2012) 

We use an extension of the van der Pol oscillator as an example of a system with multiple time 
scales to study the susceptibility of its trajectory to polynomial perturbations in the dynamics. A 
striking feature of many nonlinear, multi-parameter models is an apparently inherent insensitivity 
to large magnitude variations in certain linear combinations of parameters. This phenomenon of 
"sloppiness" is quantified by calculating the eigenvalues of the Hessian matrix of the least-squares 
cost function. These typically span many orders of magnitude. The van der Pol system is no 
exception: Perturbations in its dynamics show that most directions in parameter space weakly affect 
the limit cycle, whereas only a few directions are stiff. With this study we show that separating 
the time scales in the van der Pol system leads to a further separation of eigenvalues. Parameter 
combinations which perturb the slow manifold are stiffer and those which solely affect the jumps in 
the dynamics are sloppier. 
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INTRODUCTION 

In this manuscript, we analyze the sensitivity of a 
multiple time scales dynamical system to perturbative 
changes in its evolution laws. Rather than utilizing the 
traditional means of examining the structural stability for 
probing qualitative changes to the attractor as a response 
to perturbations, we study the structural susceptibility 
for quantifying the effects of the perturbations on the 
time series [1]. More specifically, we ask how sensitive 
is the dynamical system dz/dt = f(z) to infinitesimal 
changes of the form dz/dt = f(z) + a • g(z), for a family 
of perturbations g(z) when the parameters a ^ 0. 

This report introduces the new concept of "structural 
susceptibility" in dynamical systems, and is an outgrowth 
of our group's previous work on "sloppiness" in multipa- 
rameter systems wherein we have found that data-fitting 
in a number of nonlinear, multiparameter models is only 
sensitive to a few directions in parameter space at the 
best fit [l"!!]- The key difference between studying slop- 
piness and structural susceptibilities is that in the former, 
the parameters are intrinsic to the system, i.e., there 
are no externally introduced changes in their evolution 
laws. Nonetheless, the methodology we have developed 
for studying sloppy models is also suited for studying 
structural susceptibilities of dynamical systems. Our ap- 
proach cleanly isolates and ranks the directions in pa- 
rameter space in order of relevance to observed behavior, 
and has previously led us to suggest improvements in ex- 
perimental design [5|, extract falsifiable predictions from 
experiments [g*], and develop faster minimization algo- 
rithms [7]. Others have developed our ideas to suggest 
further improvements in experimental design (8| and pa- 
rameter estimation [9] , to quantify robustness to param- 
eter variations [10], and to set confidence regions for pre- 
dictions in multiscale models In this paper, we bring 
similar ideas together to analyze sensitivities of time se- 



ries to perturbations in dynamical systems. 

We demonstrate the utility of our approach with ap- 
plication to a dynamical s;Y^tem with two time scales — 
the van der Pol oscillator [12| which is a single parameter 
system and hence not amenable to sloppy model analy- 
sis. Instead, by choosing appropriate perturbations g(z), 
we calculate the susceptibility of its dynamics: We make 
perturbations on the attractor, and then systematically 
increase the separation of time scales in its dynamics to 
show how it can generally enhance the sloppiness in non- 
linear systems. 



MULTIPLE TIME SCALE DYNAMICS 

Multiple time scales are often found in the solutions of 
dynamical systems [13]. Broadly speaking, the defining 
criterion of these models is that the trajectory of one or 
more phase variables has an identifiable fast piece such as 
a jump or a pulse and a slow piece where the value of the 
variable doesn't change quickly [l3]. In two dimensions, 
these systems are commonly studied in the contexts of 
slow-fast vector fields written as: 



ex =X{x,y,e), 
y =Y{x,y,e) 



(1) 



where the parameter e > is small and dot indicates 
derivative with respect to time t. For 0{1) functions X 
and r, and X / 0: x = 0{l/e) and y = (9(1), so that e 
is the ratio of time scales in the system. On one extreme, 
the singular limit e = corresponds to a differential alge- 
braic system X = 0, Y = y where the solutions of X = 
comprise the "critical manifold" close to which the flow 
in phase space is slow (the "slow manifold"). Similarly, 
e = 1 corresponds to a limit where there is no separation 
of time scales, with a crossover at intermediate values of 
e. 
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Originally introduced in 1927, the van der Pol equa- 
tion, X — — x'^)x -\- X = 0, is a well- studied example 
of a second-order, nonlinear system with multiple time 
scales in its solution. Using the Lienard transformation 
y = x — x^/3 — x/jii^ and redefining time t tja, the equa- 
tion can be written as a two dimensional system 13, HI 
given by: 



/x = 100 



jii '^x 

y 



X - 



(2) 



which has the same form as ([T]) with e = The global 
attractor of this dynamical system is a structurally stable 
limit cycle with two time scales [16j . 

The van der Pol system provides a convenient way to 
separate time scales by varying ja: Small values of ja in 
the van der Pol system correspond to a small separation 
of time scales. It can be shown that the trajectory ap- 
proaches that of the harmonic oscillator as /i ^ [12]. 
At large values of /i, the system shows a separation of 
time scales which increases with increasing /i. As shown 
in figure [T] (b, c), with increasing /i, the trajectory of x 
separates into a slow part that lies 0{exp{ii~'^)) close to 
the phase space curve given by x = 0, ie. the critical 
manifold y = x — and a fast part which connects 

the two branches of the slow flow. Likewise, the separa- 
tion of time scales in y are associated with the increasing 
sharpness of the kink in its trajectory. 

The fact that with an increasing separation of time 
scales the trajectory spends an increasing amount time on 
the slow manifold and a decreasing amount of time on the 
jumps has important implications for fitting parameters 
to time series data of the van der Pol system. With 
increasing scale separation, one expects that the cost of 
fitting will be decreasingly sensitive to changes in the 
jumps of the trajectory as they get progressively shorter 
in duration. 



SLOPPINESS IN NONLINEAR FITS 

In this section, we discuss the concepts of sloppiness 
and structural susceptibility in more detail with exam- 
ples as a prelude to our calculations. For time series 
z(t,a), a least-squares fit to data di minimizes a cost 
C = I ^ii^iti^Si) — di)'^/af in the space of system pa- 
rameters which are collectively denoted as a. Our dis- 
covery of sloppiness is essentially that the eigenvalues 
of the Hessian of the cost with respect to parameters, 
= d'^C / da^da^^ at the best fit span many orders 
of magnitude. The larger and smaller eigenvalues corre- 
spond to stiffer and sloppier directions respectively. For 
concreteness, consider a time series of a multi parameter 
model, such as the one denoted by ?/(r) in figure [TJb, 
top row). The error bars schematically show the least- 
squares fit of 7/(r) and the sidebar (figure [TJa)) shows the 
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FIG. 1. (a) Eigenvalues of the Hessian matrix of the cost of 
fitting at 11 — 1, and (b, c top row) one period of time series 
x(t) and y(r), for < r < 1, are shown for ii — 1 and /x = 100 
as function of time along with schematic error bars for the 
data-fitting of the trajectory of variable y. (d) Eigenvalues for 
/i = 100 are shown on the right. As /x ^ oo, the orbit collapses 
onto the critical manifold with the trajectory spending most 
of its time on the slow manifold and vanishingly short on the 
jumps. Also shown in (b, c bottom row) is the orbit in xy 
plane (solid line) and the critical manifold (dashed line). 



eigenvalues of the corresponding Hessian matrix. Note 
the broad range of eigenvalues (~ 10^^, corresponding 
to a factor of almost a million in parameter range) — a 
feature that turns out to be typical in nonlinear fits. 

Another vivid example of sloppiness in nonlinear mod- 
els is provided by the well-established formalism behind 
the characterization of the sensitivities of initial condi- 
tions using Lyapunov exponents [17]. Consider dz^/dt = 
f(z) as a model whose parameters are the initial con- 
ditions = Zoc{^) and whose predictions are the final 
positions Zi{t) at time t. At the best fit, l-Loi(3 — {J^J)a/3 
where Jia = dzi{t)/dza{0) is the Jacobian of the sensitiv- 
ities to perturbations in the initial conditions. The Lya- 
punov exponents, which are defined to be the eigenval- 
ues of L = limt^oo l/(2t) log( J-^ J), utilize the same 
Hessian we would use in calculating the sloppy model 
eigenvalues = exp{2Un)' The roughly equal spacing 
of Lyapunov exponents naturally explains both the expo- 
nentially broad range of sloppy model exponents and the 
roughly equal spacing of log(An) for a model with initial 
conditions as parameters. 

Instead of the sensitivities with respect to the initial 
conditions or other intrinsic parameters, we focus here 
on the sensitivity of the dynamics to changes in the dy- 
namical evolution laws. Therefore, for the remainder of 
this paper we will be interested in a cost function that 
measures the square of the distance between two time 
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series for the system dz^/dt = f(z) + a • g(z) — one with 
perturbations — z(t,a 0), and the other one without, 
i.e., z(t, a=0) 



C 



2 Jo 



|z(t,a-> 0) -z(t,a = 0) 



■dt 



(3) 



with the perturbing terms gi{z) giving a power series in 
the components of z. Further in the manuscript, we wih 
use this form of the cost to compute the susceptibihty 
of the van der Pol system and show how sloppiness is 
enhanced by increasing separation of time scales in the 
van der Pol equations. This is in essence captured by 
figure [TJa & d) where we show that an increase in the 
van der Pol parameter /i from 1 to 100 produces roughly 
a million- fold increase in the spread of eigenvalues. 



SUSCEPTIBILITY OF VAN DER POL SYSTEM 

We perturb the van der Pol system in (|2]) by adding 
a series of additional terms. There is a long tradition 
in dynamical systems of studying equations of motion 
of polynomial form indeed, the theory of nor- 

mal forms suggests that general dynamical systems, even 
at bifurcations, can be generically mapped into a polyno- 
mial form by a nonlinear but smooth change of variables. 
Adding extra polynomial terms is routinely done to 'un- 
fold' the qualitative behavior near bifurcations [19]. Here 
we focus on quantitative changes far from bifurcations. 
In choosing our perturbations, we must cut off the poly- 
nomials at some order. There are two ways in which we 
specialize our general susceptibility analysis to the two 
time scale, periodic limit cycle of the van der Pol system. 
First, we choose the family of perturbations of order 37V 
as follows: 



y 



m-\-n<N 



(4) 

This choice has two noteworthy features — (a) We have 
grouped the polynomial perturbations so that, for m 7^ 
they vanish on the critical manifold, y = x — /3. That 
is, the parameters am,n with m 7^ do not significantly 
affect the dynamics on the slow manifold; we call these 
the "fast parameters" and correspondingly the ao,n are 
"slow parameters". The parameter ai^o duplicates /i to 
the same effect as changing the period, and we omit it. 
Surely, the eigenvalue spectrum of the general polynomial 
expansion, am,nX^y^^ behaves qualitatively similarly to 
the one presented here but our parametrization greatly 
simplifies the analysis of the eigenvector perturbations, 
(b) We only perturb the x equation. Our choice corre- 
sponds to a general expansion of a second-order equation, 
with the acceleration y = x written as a polynomial in 
the position y and velocity y = x. Perturbing both equa- 
tions produces similar behavior. 



Second, we modified the cost to focus on the limit cycle 
of the van der Pol system in two ways — (a) by rescaling 
all trajectories in our analysis so that they have the same 
unit period, and (b) by changing the initial condition so 
that it lies on the perturbed orbit and that the perturbed 
and unperturbed orbits are in phase with each other [2o| . 
When we correct the period T by initial conditions 
Yo t>y ^Yo^ arid do an overall rescaling of the time variable 
t tT, the cost functional for the time series of y{r) at 
each /i takes the following form: 

1 

= -J^ [y{r, a + fe, Yo + (^Yo, T ^ ST) - 

^(r,a,Yo,T)]2dr (5) 

In principle, changes in both time series, x{r) and y{r) 
could be incorporated in the cost function, but we get 
qualitatively similar results by keeping either or both 
variables. Choosing to measure changes only in y{r) cor- 
responds again to studying the second-order equation for 
y as an expansion in y and y. 




FIG. 2. Eigenvalues of Hessian matrix are shown here as a 
function of /a. The range 1 < /x < 100 corresponds to a ratio 
of time scales 1 < e < 10000. The five largest eigenvalues 
correspond to stiff directions in the parameter space: these 
directions perturb the slow manifold. The remainder affect 
the transient part of the trajectory which becomes smaller 
with an increasing separation of time scales and hence these 
directions are decreasingly relevant. 



The susceptibilities are still given by the Hessian ma- 
trix at the best fit (a=0): 



da^dap 

which can be written out more completely as: 
_ f 9y ^ dy dy^ ^ dy dT 

^ f dy ^ dy dy^ _^ 9y dT \ ^ 
\9a/3 dy^dajs dT da^s J 



(6) 
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Here, each of the two terms in the integral is to be inter- 
preted as a Jacobian matrix, a mapping from the finite 
dimensional parameter space to the infinite dimensional 
data space: 



Jtc 



%(r) , dy{T) gyp , dy{T) dT 

daa 



(7) 



djQ daa dT daa 

The sensitivity trajectories in the Jacobian, dy/daa^ 
dy/dyQ^ and dy/dT^ were computed using the open 
source SloppyCell package [2l|, [22[ . The expressions for 
the time invariant quantities, dy^/daa and dT/daa, were 
found by enforcing periodicity of the perturbed time se- 
ries denoted by y(r) = {x{r),y{r)) as follows: 

y (r = 0, a + fe, Yo + (^Yo, ^ + ST) = 
Y(r = 1, a + 5a, Yo + ^Yo. ^ + ^T), 

Taylor expansion of both sides of the previous equation 
leads to a vector equation: 

dy 



1 aa 



5a- 



from which both constants can be computed following 
the convention that the component denoting the change 
in initial condition of y{r) in ^Yo zero. Now with 

the Jacobian calculated, the Hessian at best fit is simply 

n = .Fj. 



Eigenvalues and Eigenvectors 

We computed the Hessian matrix given by the previous 
equation at multiple values of ji. The spread of eigenval- 
ues (figure [2]) increases as a function of ji confirming that 
sloppiness increases with an increasing separation of time 
scales. Not surprisingly [23|, for A/" = 4, the 14 eigenval- 
ues for /i = 1 already span 11 orders of magnitude, while 
for ji = 100, we observe that the stiffest eigenvalue is 18 
orders of magnitude larger than the smallest one — the 
spread increases by 10^ when ja increases to 100. 

Taken together with the eigenvectors shown in figure [3l 
some interesting facts come to light: Figure [2] shows that 
with increasing /i, the eigenvalues separate into two clus- 
ters of closely related decay exponents. The largest N 
eigenvalues approach constants. The other eigenvalues 
decay with power laws: two modes with exponents be- 
tween —2 and —3 and the remaining with exponents be- 
tween —5 and —6. Similarly, figure [3] shows that the 
eigenvectors also separate into two groups with increas- 
ing /i: The stiffest directions are linear combinations of 
the slow parameters whereas the sloppy directions are 
comprised of other parameters as expected. 

We can understand the effect of perturbations in pa- 
rameter combinations given by the eigenvectors (called 
eigenparameters) ek more clearly by observing their be- 
havior in the data space. The Jacobian transforma- 
tion of ([7]) projects the eigenvectors to the data space: 



5yk = J-e/c/\/A^ where Xk corresponds to the k^^ largest 
eigenvalue. Defined this way, these data space vectors, 
called eigenpredictions Q, 5yk^ are also orthonormal. Al- 
ternatively, the eigenpredictions are the left singular vec- 
tors in the singular value decomposition of the Jacobian 
(i.e. they are the columns of the unitary matrix U in 
J = U^V^ [24]) As shown in figureg]for /i = 1, 10 & 100 
(top three rows), we learn from the eigenpredictions that 
the stiff modes affect behavior both along the slow mani- 
fold and at the jumps. Moreover with increasing /i, as the 
eigenvalues associated with the stiff directions approach 
constants (figure [2]), so do the stiff eigenpredictions (fig- 
ure |4] rows 2, 3 columns (a) and (b)). The sloppy modes 
on the other hand, affect the dynamics on the jumps only. 
The maximum amplitudes of the (normalized) sloppiest 
eigenpredictions appears to increase roughly proportional 
to /i (corresponding to the jump duration of ~ /i~^). 
In the limit, these become (5- functions and derivatives 
concentrated at the jumps. Figure |4] (bottom row) also 
shows the limit cycles {eigencycles) with eigenparameter 
perturbations as phase space trajectories {x,y -\- r] 5yk) 
for small rj. 



DISCUSSION 

In this paper, we have introduced a formalism we call 
"structural susceptibility" for analyzing the quantitative 
dependence of dynamical systems to perturbations of the 
equations of motion. It is a generalization of 'unfolding' 
methods of bifurcation theory and the Lyapunov expo- 
nents governing the dependence on initial conditions, and 
exposes the ubiquitous presence of broad range of sloppy 
eigendirections in parameter space — largely unimpor- 
tant to the dynamics. We used this method to study 
the role of time scale separation in enhancing the sloppi- 
ness of the susceptibility spectrum in the particular case 
of the van der Pol oscillator. 

By extending the framework of our sloppy model anal- 
ysis to systems where changes in evolution laws are to be 
studied, our method offers a simple way to calculate the 
effects of broad classes of perturbations. By studying the 
structural susceptibility of a dynamical system with two 
time scales, the analysis presented here showed that slop- 
piness of nonlinear systems is enhanced by separation of 
time scales in the dynamics. With increasing separation 
of time scales in the van der Pol oscillator, the trajec- 
tory spends an increasing amount of time on the slow 
manifold and a vanishingly small amount of time in the 
transition region. The cost of perturbations is integrated 
over time and therefore we are unsurprised that the per- 
turbations that change the slow manifold will accrue the 
most cost and therefore manifest as stiff modes of the 
Hessian matrix. The remaining directions are sloppy as 
they only affect the behavior at the jumps or the fast 
pieces. These perturbations manifest as (5-functions and 
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FIG. 3. Hessian eigenvectors are shown for /x = 1, 10, and 100. Each colored small square shows the magnitude of an eigenvector 
component (the scale bar shown on the right). Eigenvectors for each /i are sorted so that the stiffer ones appear on the left; 
individual components are sorted so that "slow parameters" appear on the top. Note that with increasing /x, the stiff and sloppy 
eigenvectors separate by parameters: The stiff eigenvectors only have projections along the slow parameters which perturb the 
slow manifold, whereas the sloppy directions have projections along the fast parameters which mainly perturb the jumps. 



their derivatives — significantly affecting the phase-space 
trajectory, but over only the fast times asymptotically 
ignored in the least-squares cost. It remains a challenge 
to connect separation of time scales to parameter sensi- 
tivity in more complicated systems, but the analogy of 
the van der Pol system's behavior with other nonlinear 
physical systems of interest is clear. 

Many important dynamical systems have multiple time 
scales in their solutions: examples include models in neu- 
roscience (such as Hodgkin-Huxley model), systems bi- 
ology or chemical reaction systems (such as protein net- 
work models), and in engineering (such as models of com- 
bustion, lasers, locomotion, etc.). Our analysis suggests 
that any system with multiple time scales should become 
sloppier as the scales separate for the same reasons as we 
found in the van der Pol: Some parameter combinations 
will only affect the fast dynamics, and lead to sloppy 
modes. Perturbations which affect the slow dynamics 
will presumably accrue more cost and be stiff. 

More broadly, the sloppiness exposed by our struc- 
tural susceptibility analysis has clear implications for at- 
tempting to reconstruct the equations of motion from 
experimental data [lH because the true dynamics along 
any sloppy eigendirection will be relatively poorly deter- 
mined. This discovery has already influenced work on ex- 
perimental design optimization: estimating parameters 
is challenging [sl, l26|, but extracting predictions with- 
out constraining parameters is straightforward [6]. We 
further anticipate that the concept of structural suscep- 
tibility will be useful for studying systems with chaos, 
bifurcations and phase transitions; quantifying the un- 



foldings of these systems may also be useful for gaining 
a deeper understanding of the phenomena they model. 
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FIG. 4. Top three rows: Eigenpredictions Si/k for A; = 0, 3, 6, 9, 12at/i=l, lO&lOO are shown in sohd red hnes for 
stiff modes and sohd green for sloppy modes. These curves show the response of perturbations if the parameters are changed 
infinitesimaUy along the Hessian eigenvectors: A parameter change of norm e along eigendirection n will change the trajectories 
by AnC times the eigenpredicton. Dashed lines show unperturbed van der Pol solution for comparison (y scale on the right 
hand side). As the time scales separate, the amplitudes of the sloppiest eigenpredictions increase (roughly in proportion to /i) 
getting increasingly concentrated at the jumps. Bottom row shows the eigencycles for /x = 100 in solid red lines corresponding 
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whereas the sloppy modes only affect the jumps. 
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